program Small_code_Twiss
use my_own_little_code_utilities
implicit none
type(my_taylor) Z(2),r2,x,p,PHASE_ADVANCE,INVARIANT
type(my_taylor)  x2_average,x_average,x_average_xp
type(normalform) NORMAL
type(inverse_dragt_finn) idf
type(my_map) M,id,a_cs,disp,A_l,A_nl
real(dp) L,k1,fix(3),ma(2,2),h,dx_average_dj
complex(dp) zz
integer i,n,mf,j,i1,i2,k,nst
integer beta(2),twoalpha(2),gamma(2)
type(magnet) :: lattice(6*10)
type(ray) r
!   mad-x lattice
!   L : drift, L= 1.0;
!   QF   : QUADRUPOLE, L=1.0, K1=1.0;  
!   QD   : QUADRUPOLE, L=1.0, K1=-1.0; 
!   sf  : sextupole,  K2= 2.0;
!s1 : LINE=  10*(QF,sf,L,QD,sf,L);


Write(6,*) " Creating the simple lattice: 10*(QF,sf,L,QD,sf,L) "
nst=100
L=1.0_dp
h=  twopi/20/L
do i=1,10
 k=(i-1)*6
LATTICE(1+k)%NAME="QF"
 LATTICE(1+k)%L=L
 LATTICE(1+k)%BN=0.0_DP
 LATTICE(1+k)%BN(0)=h
 LATTICE(1+k)%h=h
 LATTICE(1+k)%BN(1)=1.0_DP
 LATTICE(1+k)%N=nst
LATTICE(2+k)%NAME="SF"
 LATTICE(2+k)%L=0.0_DP
 LATTICE(2+k)%BN=0.0_DP
 LATTICE(2+k)%BN(2)=2.0_DP
 LATTICE(2+k)%h=0
 LATTICE(2+k)%N=0
LATTICE(3+k)%NAME="L"
 LATTICE(3+k)%L=1.0_DP
 LATTICE(3+k)%BN=0.0_DP
 LATTICE(3+k)%h=0
 LATTICE(3+k)%N=1
LATTICE(4+k)%NAME="QD"
 LATTICE(4+k)%L=L
 LATTICE(4+k)%BN=0.0_DP
 LATTICE(4+k)%BN(1)=-1.0_DP
 LATTICE(4+k)%BN(0)=h
 LATTICE(4+k)%h=h
 LATTICE(4+k)%N=nst
LATTICE(5+k)%NAME="SF"
 LATTICE(5+k)%L=0.0_DP
 LATTICE(5+k)%BN=0.0_DP
 LATTICE(5+k)%BN(2)=2.0_DP
 LATTICE(5+k)%h=0
 LATTICE(5+k)%N=0
LATTICE(6+k)%NAME="L"
 LATTICE(6+k)%L=1.0_DP
 LATTICE(6+k)%BN=0.0_DP
 LATTICE(6+k)%h=0
 LATTICE(6+k)%N=1
enddo

delta_is_3rd_parameter=.false.
epsclean=1.e-9


MY_ORDER=4
BETA(1)=0; BETA(2)=2;  !  COEFFICIENT OF P**2 OF INVARIANT R2
GAMMA(1)=2; GAMMA(2)=0;  !  COEFFICIENT OF X**2 OF INVARIANT R2
TWOALPHA(1)=1; TWOALPHA(2)=1;  !  COEFFICIENT OF XP OF INVARIANT R2  = 2*ALPHA

X=1.0_DP.MONOMIAL.1 ;P=1.0_DP.MONOMIAL.2; R2 = X**2+P**2 ;

FIX=0.0_DP  ! FIXED POINT
CALL FIND_CLOSED_ORBIT( FIX,LATTICE,1)

 

ID=1
R%Z=FIX+ID ! INITIALIZE THE RAY AS --> RAY = FIXED POINT + IDENTITY (TAYLOR MAP)

!  COMPUTING A ONE-TURN MAP TO ORDER MY_ORDER

CALL TRACK_LATTICE( R,LATTICE,1,1)
M=R%Z    
NORMAL=M    ! NORMALIZING THE MAP
A_CS=NORMAL%A_T

CALL AVERAGE(X,A_CS,X_AVERAGE,X_AVERAGE_XP,USE_J=.TRUE.)
WRITE(6,*) " NON-LINEAR DISPERSION DEFINED AS <X> = D<X>/DJ * J +  ..."
WRITE(6,*) " TAYLOR MAP RESULTS FOR <X> = D<X>/DJ * J +  ..."
CALL PRINT(X_AVERAGE,6)








!!!!!   TRACK LATTICE FUNCTIONS !!!!

!! CANONIZE FORCES THE TRANSFORMATION INTO A SPECIAL FORM
! THE LINEAR PART IS IN COURANT-SNYDER FORM A_12=0 INCLUDING
! PARAMETERS.

  R%Z=FIX+NORMAL%A_T   ! INITIALIZE THE RAY AS --> RAY = FIXED POINT + A
  CALL CANONIZE_RAY( R%Z,PHASE_ADVANCE,A_CS,DISP,A_L,A_NL );DX_AVERAGE_DJ=0.0_DP; 
  PHASE_ADVANCE=0.0_DP
 
 INVARIANT=R2*A_CS**(-1)

DO J=1,SIZE(LATTICE)

 CALL TRACK_MAGNET( R,LATTICE(J) )

 CALL CANONIZE_RAY( R%Z,PHASE_ADVANCE,A_CS,DISP,A_L,A_NL ) 

 INVARIANT=R2*A_CS**(-1)

 DX_AVERAGE_DJ=(INVARIANT.SUB.BETA)**1.5_DP*LATTICE(J)%BN(2)/4.0_DP &
   *(-SIN(PHASE_ADVANCE*TWOPI)+SIN((PHASE_ADVANCE-NORMAL%TUNE)*TWOPI)) &
   /(1.0_DP-COS(NORMAL%TUNE*TWOPI)) + DX_AVERAGE_DJ

ENDDO

DX_AVERAGE_DJ=DX_AVERAGE_DJ*SQRT(INVARIANT.SUB.BETA)

write(6,'(a10,E20.13,a20)') '  d<x>/dJ ', dx_average_dj, " < ---- analytical  "  
!  Non-Linear dispersion defined as <x> = d<x>/dJ * J +  ...



end program Small_code_Twiss

subroutine test_sum
use my_own_little_code_utilities
implicit none
complex(dp) sum
real(dp) nu,thetap,nup
integer i,n
write(6,*) "n "
read(5,*) n
!n=100
nu=0.1745634d0
101 write(6,*) "nup "
read(5,*) nup
if(nup<0) goto 100

thetap=twopi*nup
sum=0

do i=-n,n

    sum=sum-exp(-i_*i*thetap)/(i-nu)/twopi/i_
    
enddo

write(6,*) (sum)
write(6,*) (exp(-i_*thetap*nu)/(1.d0-exp(-i_*twopi*nu)))
write(6,*) abs(sum)
write(6,*) abs(exp(-i_*thetap*nu)/(1.d0-exp(-i_*twopi*nu)))

goto 101
100 return
end subroutine test_sum